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AN  HP- ADAPTIVE  METHOD  IN  SPACE  AND  TIME 


FOR  PARABOLIC  SYSTEMS* 

JOSEPH  E.  FLAHERTYt  AND  PETER  K.  MOORE* 


Abstract.  We  describe  an  adaptive  method-of-lines  ^-refinement  algorithm  in 
space  and  time  for  one-dimensional  vector  systems  of  parabolic  partial  differential 
equations.  Solutions  are  calculated  using  Galerkin’s  method  with  a  piecewise- 
polynomial  hierarchical  basis  in  space  and  singly-implicit  Runge-Kutta  (SIRK) 
methods  in  time.  A  posteriori  estimates  of  the  local  spatial  and  temporal  discretization 
error  are  used  with  a  priori  error  estimates  to  control  spatial  and  temporal  enrichment. 
Computational  results  are  used  to  compare  and  verify  the  utility  of  several  variants  of 
the  basic  /z/t-refinement  procedure. 

Key  words,  adaptive  refinement,  finite-element  methods,  a  posteriori  error  estima¬ 
tion. 


MSC  subject  classifications.  65M20,  65M50,  65M60 

1.  Introduction.  Basic  adaptive  strategies  for  the  solution  of  parabolic  partial 
differential  equations  (PDEs)  include  method-of-lines  (MOL)  approaches  with  spatial 
mesh  refinement  or  coarsening  (/z-refinement)  [1,  3,  5,  13,  18],  mesh  motion  (r- 
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refinement)  [1,  16],  and  order  variation  (p-refinement)  [3,  14,  18].  Solution  “enrich¬ 
ment  indicators,”  which  are  frequently  estimates  of  local  spatial  discretization  errors 
[1-3,  5,  18,  21],  are  obtained  from  preliminary  solutions  and  used  to  identify  portions 
of  the  domain  in  need  of  additional  resolution.  Some  combination  of  the  three  basic 
enrichment  strategies  is  used  to  alter  the  approximation  and  recursively  calculate 
improved  solutions  until  specified  accuracy  criteria  have  been  satisfied.  The  focus  to 
date  has  been  on  h-  and  r-refinement  [1-3,  5,  13,  16,  18]  with  little  attention  devoted 
to  p-  and  /zp-refinement  [3,  19].  Since  the  latter  schemes  have  been  remarkably  suc¬ 
cessful  for  elliptic  problems  it  is  natural  to  examine  their  utility  when  solving  parabolic 
systems.  With  few  exceptions  [17],  nothing  has  been  done  to  coordinate  spatial  and 
temporal  enrichment  and  we  address  a  combined  space-time  adaptive  /zp-refinement 


strategy  herein. 

We  consider  vector  systems  of  M  parabolic  equations  of  the  form 

u,  +  f(x,t,u,ux)  =  D(x,t,u,ux)x,  reQ  =  ( c,d ),  t  >  0,  (1.1a) 

u(x,0)  =  u0(x),  x  £  Cl,  (1.1b) 

Ui(x ,t)  =  ct(t)  x  £  dClf,  (1.1c) 

Di(x,t,u,ux)  =  x£dCl[f,  i  =  1,  2,  ••*,  M ,  t  >  0.  (l.ld) 


The  boundary  dCl  =  dCl f  +  dCl?  is  divided  component- wise,  i  =  1,  2,  M ,  into  sets 
where  essential  ( E)  and  natural  (AO  data  is  applied.  Additional  restrictions  must  be 
placed  on  the  functions  f  and  D  to  ensure  that  (1.1)  is  a  well-posed  parabolic  problem 
with  a  locally-isolated  solution. 

We  solve  (1.1)  using  a  finite  element  Galerkin  MOL  technique  with  a  piecewise 
polynomial  hierarchical  spatial  basis  of  degree  p  ^  1  and  a  singly  implicit  Runge- 
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Kutta  (SIRK)  method  in  time  (§2).  SIRKs  are  well-suited  to  adaptive  computation 
which  requires  frequent  restarts  of  the  temporal  integration  due  to  spatial  enrichment,  a 
consideration  that  would  be  even  more  important  were  local  temporal  refinement  [18] 
used  instead  of  a  MOL.  SIRKs  additionally  have  high  stage  order  that  eliminates 
order  reduction  [11,  21]  and  error  estimates  of  each  stage  that,  with  stability  considera¬ 
tions,  permit  acceptance  of  solutions  at  intermediate  stages  whenever  the  final  solution 
lacks  the  requisite  accuracy  [20]. 

A  posteriori  temporal  (§3.1)  and  spatial  (§3.2)  error  estimates  are  used  to  guide 
the  adaptive  /zp-refinement  algorithm  as  described  in  §4.  A  temporal  order  and  step 
selection  strategy  (§4.1)  is  used  to  develop  a  base  adaptive  /zp-refinement  strategy  and 
several  variants  (§4.2).  Four  examples  are  presented  in  §5.  A  comparison  of  the  spa¬ 
tial  /zp-spatial  refinement  strategies  of  §4,  an  /2-refinement  strategy  with  various  fixed 
orders  (§4.2),  and  a  p  -refinement  strategy  with  uniform  spacing  is  made  using  a  linear 
example.  Burgers’  equation  is  used  to  compare  strategies  for  determining  initial 
guesses  for  Newton’s  method  and  for  comparing  spatial  /zp-refinement  strategies.  A 
comparison  of  spatial  /zp-refinement  strategies  is  also  made  using  a  Brusselator  with 
some  slight  diffusion.  The  Brusselator  [15]  and  a  shear  band  model  [12]  demonstrate 
the  effectiveness  of  our  method  on  realistic  nonlinear  systems.  A  brief  discussion  of 
our  results  is  presented  in  §6. 

2.  Discretization.  The  Galerkin  form  of  (1.1)  consists  of  determining 
u (x,t)£  He(Q) x (t  >  0)  such  that 

(u,  ,v)  +  (f,v)  +  (D.v* )  =  Dr  v  I  ,  for  all  v  £  H$ ,  t  >  0,  (2.1a) 

x  e  an 

(u,v)  =  (u0,v),  for  all  v  G  H$ ,  t  =  0,  (2.1b) 

where 
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(u,v)  =  juT\dx.  (2.1c) 

As  usual,  the  Sobolev  space  H1(Q)  consists  of  functions  having  square  integrate  first 
derivatives  and  the  subscripts  E  and  0  further  restrict  functions  to  satisfy  (1.1c)  and  a 
trivial  version  of  (1.1c),  respectively. 


Introduce  a  partition 


A q  " —  { c  —  x  q  <  x  ^  <  ■■■  K Xft  —  d }  (2.2) 

of  Q  into  N  subintervals  and  approximate  Hl  by  a  finite-dimensional  subspace  SAa 
where  SAq  consists  of  piecewise  polynomials  whose  restriction  to  (xi_1^ci)  are  polyno¬ 
mials  of  degree  pt  ^  1,  i  =  1,  2,  •••,  N.  Consider  a  finite  element  approximation 

U(x,/)G  SAq  of  u (x,t)£  He  having  the  form 

U(x,0  =  LfWWnC*)  +  £  £U/*(')<M*)  (2.3a) 

j=0  i=lk=2 


where  here,  and  throughout  this  paper,  Galerkin  coordinates  are  identified  by  an  over¬ 
bar.  The  functions 


M*)  =  i 


(x  -  *,_!)/(*,•  -  xt_{  <,x  <xt 

(xi+ 1  -  x)/(xi+1  -  x*),  Xj  (  £  X  <  xl+l,  i  =  0,  1,  •••,  N,  k  =  1, 
0,  otherwise 


(2.3b) 


and 


x 


<M*)  = 


I  J  ^/,*-i(v)^v. 


1  ^  X  <  Xi 

otherwise 


i  =  1,  2,  •••,  N,  k  =  2,  3,  •  Pi, 


(2.3c) 


where  /^(x)  is  the  k th  -degree  Legendre  polynomial  scaled  to  the  subinterval 
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comprise  a  hierarchical  basis  for  SA°  [25]. 

In  a  similar  manner,  we  approximate  v  by  VG  SAq;  thus,  replacing  u  and  v  in 
(2.1)  by  U  and  V,  we  determine  U(x,0  as  the  solution  of  the  ordinary  differential  sys¬ 
tem 

(U, , V)  +  (f, V)  +  (D,VX )  =  DrV  I  ,  for  all  V G  S 0An ,  t  >  0,  (2.4a) 

x  e  do 

(U,V)  =  (u0,V),  for  all  V  G  S  0An ,  t=  0.  (2.4b) 

Explicit  utilization  of  (2.3)  reveals  that  (2.4)  has  the  form  [20] 

MU"  =  F(U),  t  >  0,  MU(0)  =  h(uo),  (2.5a,b) 

where  M,  F,  and  h  are  the  mass  matrix,  load  vector,  and  the  initial  load  vector,  respec- 

_  N 

tively;  U(r)  is  a  vector  of  Galerkin  coordinates  of  length  n  =  M(N  +  1  +  £(p*  -  1)) 

i=l 

on  Sgn;  and  ( )'  denotes  total  differentiation. 

Implicit  numerical  methods  are  generally  preferred  for  the  temporal  integration  of 
systems  like  (2.5)  that  arise  from  parabolic  problems.  Stiffness  arises  in  a  MOL  for¬ 
mulation  due  to  the  spatial  discretization,  so  backward  difference  formulas  (BDFs)  are 
a  common  solution  technique  [1-3,  5,  17].  SIRKs,  introduced  by  Butcher  [9]  and  Nor- 
sett  [22]  and  extended  by  Burrage  [7],  offer  A  -stability  in  combination  with  high-stage 
order  and  Jacobians  of  equal  dimension  to  those  of  multistep  methods;  hence,  they 
may  provide  an  attractive  alternative  to  BDF  methods.  We  depart  from  Butcher’s  [9] 
SIRK  procedure  by  approximating  U(r)  rather  than  U'(r)  for  the  time  step  [a, a  +  h] 
by  an  expression  of  the  form 

U(/)«W(0=£W,v,((),  v,(0= 

1=0  m= 0  " 


(2.6a,b) 
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where  the  coefficients  alm,  l,  m  =0,  1,  s,  axe  given  by  Moore  and  Flaherty  [20], 
W,  =W(a+Clh),  1  =  1,2,  •••,  s,  and  W0  =  W(a).  Spatially-dependent  functions 
W(*,0  and  W,(x),  /=0,  1,  •••,  s,  may  be  defined  by  using  W(r)  and  W,, 
/  =0,  1,  •••,  s,  respectively,  in  conjunction  with  (2.3a-c). 


The  parameters  cz  are  selected  as  A^,,  /  =  1,  2,  s,  where  ^  is  the  Ith  root  of 
Ls(t),  the  Laguerre  polynomial  of  degree  5.  As  described  by  Burrage  [7],  choosing  X 
as  1/^*,  1  ^  1*  <,  s,  can  produce  a  method  with  favorable  stability  and  discretization 
error  properties.  Using  properties  of  Laguerre  polynomials,  Butcher  [9]  introduced  the 
transformation 


[W1,W2,-,W,]r  =  T[W1,W2,— ,W5]^ 


where  T  is  the  tensor  product  of  T, 


(2.7a) 


Tu  i,  j  =  1,  2,  s,  (2.7b) 

with  an  n  x  n  identity  matrix  that  reduces  the  linearized  version  of  (2.5a,  6)  to  the 
form  [20] 


M-MFjj  0 

•  •  0 

< 

M  M-A/iFd 

•  •  0 

aw2 

h2 

M  M 

•  •  M  -  AAFp 

AW, 

Hi 

h2 

MW0 

MW0 

M  0  •  •  ■  0 

M  M  0 

'  »-H  CN 

A 

Fi 

f2 

H, 

MW0 

+ 

M  M  •  •  •  M 
_  _ 

w,_ 

-  Xh 

F 

x  s 

Thus,  the  Newton  corrections  AWZ  to  the  transformed  variables  Wz  ,/  = 


(2.7c) 


(2.7d) 


1,  2,  •••,  s. 
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are  determined  by  a  simple  block  backward  substitution  involving  s  linear  systems. 
Only  one  Jacobian  is  involved  in  all  Runge-Kutta  stages  and  inner  products  involving 
the  lower  diagonal  entries  in  (2.7c, d)  may  be  accumulated  stage  by  stage.  Addition¬ 
ally,  the  terms  MW;,  /  =  1,  2,  •••,  s,  in  (2.7 d),  need  not  be  evaluated  after  the  initial 
Newton  iteration  since  corrections  of  MAW,,  /  =  1,  2,  •••,  s,  are  generated  by  (2.7c). 
Inverting  (2.7b)  to  obtain  W,,  /  =  1,  2,  s ,  is  straightforward  [9],  By  using  the  for¬ 
mulation  (2.6-7)  instead  of  Butcher’s  [9]  original  approach,  we  avoid  solving  an  addi¬ 
tional  linear  system  by  determining  W  directly  rather  than  obtaining  it  from  F. 

With  X  =  1  fer,  1  £  /*  £  5,  we  have  cr  =  1;  hence,  W(a+/z)  =  Wr  is  regarded 
as  the  final-stage  solution  that  is  propagated  to  the  next  time  step  if  accuracy  condi¬ 
tions  are  satisfied.  All  other  solutions  W/t  /  =  1,  2,  — ,  /*  -  1,  /*  +  i,  ^  furnish 

intermediate  results.  Moore  and  Flaherty’s  [20]  stability  analysis  permits  the  accep¬ 
tance  of  solutions  at  some  intermediate  stages  for  /  <  /*.  This  is  particularly  useful 
when  the  discretization  error  of  W,*  is  larger  than  prescribed  but  that  of  W,,  /  <  /*,  is 
acceptable.  Those  stages  that  satisfy  Moore  and  Flaherty’s  [20]  stability  conditions 
and,  hence,  may  be  propagated  forward  are  called  “acceptable  stages.”  Acceptable 
stages  with  /  <  /  and,  hence,  c,  <  1  are  refenred  to  as  “partial  steps”  while  “full 
steps  are  taken  when  1=1.  For  reference,  we  include  the  set  of  acceptable  stages, 
Fy,  s  =  1,  2,  •••,  8,  in  Table  2.1  [20],  For  the  adaptive  algorithm  of  §4.2  it  will  also 
be  useful  to  define  the  set  of  “adaptive  stages”  *  =  1,  2,  -,  8,  as  presented  in 
Table  2.2. 


s 

i 

2 

3 

4 

5 

6 

7 

8 

rs 

{i} 

U,2] 

{1,2} 

{1,2} 

{2,3} 

{2,3} 

{2,3,4} 

{2,3,4} 

Table  2.1.  The  set  of  potentially  acceptable  stages  for  a  s  -stage  SIRK 
s  =  1,  2,  •••.  8. 
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s 

1 

2 

3 

4 

5 

6 

7 

8 

% 

{1} 

{1,2} 

{1,2,3} 

{1,2,3} 

{2,3,4} 

Table  2.2.  The  set  x¥s  of  adaptive  stages  for  a  s  -stage  SIRK,  s  =  1,  2,  8. 


3.  Error  Estimation.  The  adaptive  algorithms  of  §4  utilize  estimates  of  the  total 
error  e  =  u  -  W,  the  spatial  error  ex  =  U  -  W,  and  the  temporal  error  e?  =  u  -  U. 
Procedures  for  estimating  the  total  error  [2,  21]  are  reviewed  briefly  for  completeness. 
All  estimates  are  calculated  using  the  root-mean-square  norm  [6] 


VI  M 


Ik 


(3.1) 


(atol1  +  rfo/'||WI||1)2 

where  el  is  the  ith  component  of  e,  and  axol1  and  rtoV ,  i  =  1,  2,  M,  are 
prescribed  absolute  and  relative  error  tolerances,  respectively.  Temporal  error  estima¬ 
tion  is  the  subject  of  §3.1  and  spatial  error  estimation  follows  in  §3.2. 


3.1.  Temporal  Error  Estimation.  Temporal  error  estimates,  which  control  step 
and  order  selection,  can  be  obtained  for  all  stages  [8,  20]  by  embedding.  The  addi¬ 
tional  stage 


Wi+1  =  W0  -  £[c,+1L't(cJ+1M)/*  +  ^-!(c,+1A)]  [£  l/A,Wf  -  WJ  +  hX¥(Ws+l), 

k= 1  i=l 

(3.2) 


furnishes  higher-order  temporal  solutions  that  may  be  subtracted  from  existing  solu¬ 
tions  to  obtain  error  estimates 


le'C^+Q/OlU  « IIE/(-)IU  »  II4+i(qA)(Wj+i0  -  W0(-))IU, 


/  =  1,  2. 


(3.3) 
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of  all  stages  [20]. 

Estimates  of  the  scaled  derivatives  l|/ri+*8/+*U(-,a  +  h)\\rms,  k  =-l,  0,  1,  2,  are 
needed  for  the  order  selection  strategy  described  in  §4.1.  Estimates  of 
/z*-15/_1U(-,<2  +  h)  and  hsd?\J(-,a  +  h)  are  obtained  by  writing  (2.6)  in  the  form  [20] 

W(0  =  Wo Lsm)  +  t  WJL^Ce/X)  -  Lsm)l  6  =  ~~  (3.4a,b) 

i=l  « 

and  differentiating  s  -  1  and  s  times,  respectively. 


Results  of  Burrage  et  al.  [8]  imply  that 

E ,'•«)  =  -h  it  (bj  -  Fy)F(Wj)  -  6s+lF(W5„)]  (3.5a) 

>1 

where  bj,  j  =  1,  2,  •••,  s,  and  bj,  j  =  1,  2,  •••,  s  +  1,  satisfy  the  consistency  condi¬ 
tions 


2  bjcj  = 


7= i 


/  +  r 


/  =  o,  1,  •••,  j-1. 


$+1  _ 

£  Vl 

7=1 


/  +  1 


,  /  =  0,  1,  5. 


(3.5b,c) 


Using  Taylor’s  series  and  (3.5b,c),  it  follows  that 


ME  (bj  ~  bj)z(a  +  cjh )  -  6J+1z(a  +  c4+1/i)] 

7=1 

=  z<-)(o)**«2r[£/,;<J-7lT]  (3.6) 

for  any  function  z(i)£CI+1[fl^i].  Extending  (3.6)  to  vector  systems  with  z  =  U' 
yields  [8] 

Hs+lVj?+1\a  +h)  =  — - — - +  O  (Hs+2)  (3.7) 

tbjc]-ll(s  +  1) 

7=1 

where  H  =  Xh . 

The  final  scaled  derivative  Hs+2\}j*s+2\a  +  h)  is  estimated  by  taking  two  full 
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steps  on  the  same  spatial  mesh  with  the  same  temporal  order  and  time  step  and  using 
(3.7)  to  show  that  [8] 


where 


Hs+2U[?+2\a 


+  h)  = 


\s+2VEj*(a  +h)sl 

t  bjc )  -  1  !{s  +  1) 
j= i 


+  0  ( Hs+ 3) 


(3.8a) 


VE/*(a  +  h)  —  E  I"  (a  +  h)  —  E  (3.8b) 

3.2  Spatial  Error  Estimation.  Estimates  of  the  spatial  and  total  discretization 
errors  are  obtained  by  computing  two  additional  solutions 

Yk(x,t)  =  W(x,t)  +  Ek(x,t),  Ek(x,t)  =  E‘(x,t)  +  ZWJ(x,t),  k  =  0,1, 

/= o 

(3.9a,b) 

where  E*  and  E^^  are,  respectively,  total  and  spatial  error  estimates  of  each  solution 
Y*,  k  =0,  1.  The  temporal  error  estimate  Ef ,  obtained  as  described  in  §3.1,  is  used 
for  all  acceptable  SIRK  stages,  otherwise  it  is  set  to  zero.  The  complexity  of  obtaining 
spatial  error  estimates  may  be  reduced  by  using  nodal  superconvergence  [2,  21],  Thus, 
errors  on  An  are  neglected  and  spatial  errors  are  approximated  by  local  p-refinement 
with 


N 


Ex'k{x,t)  =  £  #•*(*) 4^+jOO,  k  =  0,  1, 

i=l 


(3.10a) 


by  replacing  u  in  (2.1)  by  Y*,  k  =0,  1,  to  obtain  the  local  Galerkin  problems 
(Y*V)t  +  (f(x,t,Yk,Yk)y)i  +  (D(x  ,t  ,Yk  ,Y  k)i  =  0,  for  all  V  =  s^+k+Ux^\ 

t  £  (a, a  +  csh),  i  =  l,  2,  -,N,  (3.10b) 


Yk(x,a)  =  W(x,a),  k  =0,1. 


(3.10c) 


The  problems  (3.10b,c)  are  solved  for  Y*,  k  =  0,  1,  using  the  same  s  -stage  SIRK  that 
was  used  to  obtain  the  finite  element  solution  W.  The  local  spaces  ^+k+1'(,Xt  x‘-^ 
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consists  of  polynomials  having  maximal  degree  pt+k+l,  k-  0,1,  on  ,*,•), 
i  =  1,  2,  -,  N .  The  local  L2  inner  product  (u,v);  is  defined  according  to  (2.1c)  with 

£2  replaced  by  (xi_1  ,%,■). 

The  total  euor  ||e(-,a  +  Cjh)\\rms,  is  approximated  by  ||E°(-,a  +  Cjh)\\rm,  j  £  Ts; 
however,  the  spatial  order  selection  strategy  (cf.  §4.2)  requires  estimates  of 
||E^(-,a  +  Cjh)\\rms  i,  k  =  -2,  -1,  0,  1,  j  £  'Ey,  i  =  1,  2,  •••,  N.  The  estimates  E*’0 
and  E*’1  are  computed  using  (3.9,10).  The  errors  E*’-1  and  Ex~2  of  solutions  one  and 
two  degrees  lower,  respectively,  than  the  current  approximation  on  each  element,  are 
obtained  naturally  from  the  pfh  and  pt  -  Ist  terms  of  the  hierarchical  series  (cf.,  e.g. 
(2.3)).  This  approach  differs  from  the  order  extrapolation  techniques  for  elliptic  prob¬ 
lems  of  Szabo  [26]  and  Zienkiewicz  et  al.  [29].  Unlike  these  procedures,  availability 
of  Y1  provides  us  with  a  rational  basis  to  change  orders  from  piecewise  linear  to 
piecewise  quadratic  solutions. 

4,  Adaptive  Strategies.  A  top-level  pseudo-C  description  of  an  adaptive  MOL 
hp-refinement  algorithm  for  solving  (2.1)  on  Q.x(0,T]  is  presented  in  Figure  4.1.  Input 
to  this  procedure,  called  MOLhp,  consists  of  T  and  absolute  and  relative  error  toler¬ 
ance  vectors  atol  and  rtol,  respectively.  The  core  of  MOLhp  is  the  integration  of  (2.1) 
for  a  single  time  step  of  duration  h.  The  temporal  (§4.1)  and  spatial  (§4.2)  hp- 
refinement  strategies  are  based  on  the  order  and  step  selection  strategy  used  within  the 
BDF  code  DASSL  [6].  With  both  temporal  and  spatial  enrichment,  the  order  is  chosen 
prior  to  selecting  the  step  size. 

4.1  Adaptive  Strategies  in  Time.  Current  techniques  for  selecting  the  temporal 
order  and  step  size  (cf.,  e.g.,  Brenan  et  al.  [6])  on  a  subinterval  {a  ,a  +h  ]  utilize  esti¬ 
mates  of  the  four  scaled  derivatives  ||/ji+*d/+*U(-,a)||/ms ,  k  =  -1,  0,  1,  2,  where  s  is 
the  order  of  the  previous  step  (cf.  §3.1).  The  order  is  decreased  or  increased  by  one  if 
the  sequence  of  scaled  derivatives  at  different  orders  is,  respectively,  decreasing  or 
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MOLhp(ato\,  rtol,  T) 

{ 

a  -  0; 

Generate  an  initial  mesh  Aa  and  select  an  initial  time  step  with  s  =  1; 
while  (a  +h  <,T) 

{ 

Solve  (2.5)  on  (a,a+h]xAn  using  a  s -stage  SIRK  to  obtain  the  ap¬ 
proximate  solution  W(x,a+h),  the  total  error  estimates  E °(x,a  +  eft), 
l  €  Ts,  and  the  spatial  enror  estimates  Ex,k(x,a  +  cth), 
k  =-2, -1,  0,  1,1  £ 

Use  E X*(x,a  +  cth),  k  =  -2,  -1,  0,  1,  /  G  'P,,  to  obtain  the  elemental 
spatial  error  indicators  If*,  i  -  1,  2,  •••,  N,  k  =  -2,  -1,  0,  1; 

Determine  /max  =  max  /  such  that  ||E/^J|rmj  <  1; 

If  ('max  >  0) 
a  += 

Calculate  a  new  time  step  h  and  SIRK  order  s ; 

if  (I max  <  I* ) 

Generate  a  new  grid  An  based  on  the  refinement  indicators; 
else 

Generate  a  new  grid  Aa  if  significant  refinement  or  coarsening  is 
needed; 

} 

} 

Figure  4.1.  Pseudo-C  description  of  the  adaptive  MOL  algorithm  with  local 
hp-refinement  for  solving  (2.1)  on  Ox (0,7]. 


increasing  [6],  Once  the  order  has  been  selected,  the  appropriate  scaled  derivative  is 
used  with  the  local  discretization  error  formula  to  select  the  time  step  so  that  the  tem¬ 
poral  error  estimate  is  less  than  1/10  of  the  total  error  estimate  [24],  A  factor  of  1/10 
produced  much  more  reliable  results  (cf.  Example  5.4)  than  our  previous  factor  of  1/3 
[19].  If  the  temporal  error  estimate  is  too  large,  spatial  error  estimates  may  be  unreli¬ 
able  due  to  inaccurate  integration  of  (3.10).  Thus,  the  time  step  is  reduced  without 
spatial  adaptivity  when  the  temporal  error  estimate  is  more  than  20%  of  the  total  error 
estimate  and  the  total  error  estimate  is  more  than  0.8 


-  13  - 


Interpolation  provides  estimates  of  the  solution  at  times  other  than  a+  cth, 
l  =  1,  2,  •••,  s,  which  may,  e.g.,  be  used  to  compute  initial  guesses  for  Newton’s  itera¬ 
tion.  The  interpolation  strategy  used  in  the  SIRK  code  STRIDE  [8]  satisfies  (3.4)  [20] 
and  is  denoted  as  WgBC  (r).  High-order  SIRK  formulas  require  function  evaluations 
at  times  that  are  well  in  advance  of  t  =  a  +  h  for  the  higher  stages.  Thus,  divergence 
of  Newton’s  iteration  (cf.  (2.7))  can  be  expected  without  close  initial  guesses.  This 
difficulty  prompted  Moore  and  Flaherty  [20]  to  consider  the  reduced-order  formula 

W MF(a  +  6h)  =  £  W (4.1) 

*=i 

Unfortunately,  neither  W BBC(0  nor  W MF(t)  produced  reliable  solution  estimates  at 
advanced  SIRK  stages.  Solutions  W (t)  at  t  =a  +  cth,  l  =  1,  2,  •••,  s,  appear  to  con¬ 
verge  more  rapidly  than  elsewhere  on  t  >  a .  This  suggests  a  modified  Butcher  (MB) 
strategy  which  generates  initial  guesses  by  (3.4)  if  0  £  1  and  by  the  solution  at  the 
closest  SIRK  stage  of  previous  step  if  0  >  1.  Results  of  Moore  and  Flaherty  [20]  sug¬ 
gest  that  Wssc  is  more  accurate  than  WMF  when  0  £  1  and  that  the  opposite  is  true 
when  0  >  1  and  we  call  this  the  combination  (CO)  strategy.  We  use  the  initial  data 
W0  to  provide  Newton  guesses  for  each  stage  as  a  benchmark  (IC)  strategy.  A  com¬ 
parison  of  these  strategies  appears  in  Example  5.2. 

4.2  Adaptive  Strategies  in  Space.  The  solution  of  (2.7),  the  error  estimates 
||E°(-,a  +  ct  h  )||/.OTS ,  /€  Ts,  and  the  elemental  spatial  error  estimates 
||E x'k(-,a  +  Cih)\\rmsj,  l  £  'Ey,  i  =  1,  2,  N,  k  =  — 2, -1,  0,  1  are  computed  for 
each  time  step.  The  largest  “acceptable  stage,”  /maxe  Ty,  for  which 
l|E°(-,<3  +  clmah)\\rms  <.  1  is  determined  and  the  time  step  is  advanced  accordingly;  oth¬ 
erwise  the  step  is  rejected  (cf.  Figure  4.1).  A  new  grid  is  always  generated  when  a 
step  is  rejected  or  when  a  partial  step  is  taken.  A  new  grid  may  also  be  generated 
with  a  full  step  when  significant  refinement  and/or  coarsening  is  indicated.  This  must 
be  done  carefully  since  frequent  regridding  requires  additional  assemblies  and 
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factorizations  of  functions  and  Jacobians.  It  also  reduces  the  possibility  of  increasing 
the  temporal  order  since  such  increases  require  two  successive  accepted  steps  with  the 
same  grid.  However,  only  regridding  when  a  step  is  rejected  or  a  partial  step  is  taken 
results  in  too  many  rejected  steps  and  delays  mesh  coarsening  (cf.  Examples  5. 1-5.3). 
Our  compromise  is  to  regrid  on  an  accepted  step  when  more  than  50%  of  the  elements 
should  be  coarsened  or  more  than  5%  should  be  refined. 

The  spatial  error  estimates  ||Ex’*(-,a  +  l  €  'FJ,  i  =  1,  2,  ••*,  N, 

k  =  -2,  -1,  0,  1,  are  used  to  generate  elemental  refinement  indicators  If*, 
i  =  1,  2,  •••,  N,  k  =  -2,  -1,  0,  1,  as  described  below.  Once  new  orders  pt  have  been 

selected  for  each  element,  the  mesh  is  modified  so  that  if^1  Pi 

~  1/V/V  ,  i  =  1,  2,  •••,  N ,  and  is,  thus,  equidistributed  over  the  mesh.  Specifically,  each 
element  is  refined  by  2J  where  /  is  the  smallest  integer  larger  than 

llPilogii^NI-'^  Pl)  for  refinement  and  /  = -1  for  coarsening.  The  mesh  is 
represented  as  a  binary  tree  with  finer  elements  created  by  refinement  regarded  as 
offspring  of  coarser  ones.  Coarsening  of  two  adjacent  elements  at  the  same  refinement 
level  occurs  only  when  they  have  the  same  parent  and  only  when  the  local  spatial  indi¬ 
cator  associated  with  each  element  i  is  less  than  l/(5x2/?i-1ViV  ). 

The  order  selection  strategy  on  element  i  is  based  on  the  size  of  If’0.  If 
If'0  >  1/V77,  the  order  is  changed  if  the  sequence  If*,  k  =  -2,  -1,  0,  1  is  increasing 
or  decreasing.  If  If’0  <  l/(5xV/V" 2Pi)  and  if  If’~l  <  l/(5xV/V  2Pi)  the  order  is  too  high 
and  may  be  reduced  by  one,  i.e,  p-t  =  pt  -  1.  Finally,  if  the  error  on  element  i  is 
“satisfactory”  but  either  If*  <  ll(10xJN  2Pi)  or  If'~l  <  0.5/V/V  the  order  pt  is, 
respectively,  increased  or  decreased  by  one.  In  all  other  cases,  the  order  is  unchanged. 
An  increase  in  order  should  result  in  a  coarser  mesh  which  should  reduce  the  dimen¬ 
sion  of  the  discrete  system. 
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To  ensure  a  “smooth”  mesh  gradation  a  1 -neighbor  rule  [4]  is  invoked  with 
respect  to  step  size  and  order;  thus  neighboring  elements  must  differ  in  size  by  at  most 
a  factor  of  two  and  in  order  by  at  most  one.  The  order  smoothing  was  essential  in 
solving  Example  5.4. 

The  above  spatial  refinement  strategy  with  If*  =  ||Ex^(-,a  +  cph)^rms 
k  =  -2,  -1,  0,  1,  is  referred  to  as  the  base  strategy  and  is  identified  by  the  letter  B. 
Several  variations  of  Strategy  B  follow  and  these  are  compared  in  §5. 

NG  No  Grid-Prediction  Strategy.  The  grid  is  refined  only  when  the  prescribed  error 
tolerance  is  exceeded. 

FP  Full  Prediction  Strategy.  The  idea  of  predicting  where  spatial  refinement  will  be 
needed  in  subsequent  time  steps  was  introduced  by  Bietermann  and  Babuska  [5] 
who  utilized  pattern  recognition  techniques  and  the  extrapolation  of  data  at  previ¬ 
ous  times  to  predict  future  meshes.  With  SIRKs,  solution  information  is  available 
in  advance  of  the  accepted  time  step.  In  the  full  prediction  strategy,  we  utilize 
spatial  errors  of  other  SIRK  stages  to  select  a  grid  for  a  subsequent  time  step.  In 
particular,  we  set  If*  =  ^  max  ^||Ex,*(-,a  +  clh)\\rms4  if  a  step  is  rejected  or  a 

partial  step  is  taken  and  If*  =  ^  max  |jE*^(-,a  +clh)\\rms4  if  a  full  step  is 

taken,  k  =  -2,  — 1,  0,  1. 

NP  No  Partial-Step  Strategy.  Butcher  [10]  indicated  that  partial-step  capabilities 
would  be  dropped  from  future  versions  of  the  SIRK  code  STRIDE  because  of  its 
limited  benefits  and  we  seek  to  determine  whether  or  not  this  applies  in  a  MOL 
setting.  Hence,  Strategy  NP  uses  Strategy  FP  but  partial  steps  are  not  allowed. 

SX  Szabo’s  Extrapolation  Strategy.  Szabo  [26]  uses  extrapolation  of  the  spatial 
errors  of  solutions  of  local  degrees  pit  pt  -l,  and  pt  -  2  to  construct  an  estimate 
of  the  error  for  an  approximation  of  degree  pt  +  1.  This  is  straightforward  if 
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pt  >  2.  When  pt  =  2,  local  solutions  having  degrees  2  and  1  are  linearly  extrapo¬ 
lated.  When  pt  =  1,  we  could  either  (i)  obtain  a  “higher-order”  solution  by 
adding  the  spatial  eiTor  estimate  (cf.  §3.2)  to  the  linear  solution  on  element  i  or 
(ii)  construct  a  “lower-order”  piecewise  constant  (pf  =  0)  solution  using  the 
value  of  the  piecewise -linear-solution  at  the  center  ,  of  subinterval  i.  Linear  extra¬ 
polation  is  used  in  both  cases.  Results  of  both  strategies  were  very  similar  and 
computations  in  §5  were  obtained  using  the  latter  technique.  We  use  SX  with 
full  prediction  since  this  was  significantly  more  efficient  than  strategies  with  no 
prediction. 

These  /zp-refinement  strategies  are  compared  with  several  /z-refinement  strategies 
and  a  p-refinement  strategy.  The  /z-refinement  strategies  use  the  Strategy  B  with  p 
fixed  at  2,  3,  4,  and  12.  We  refer  to  these  strategies  by  their  value  of  p.  The  p- 
refinement  strategy  uses  a  fixed  uniform  grid  with  p  varying  according  to  Strategy  B 
and  the  proviso  that  p  is  increased  by  one  when  mesh  refinement  is  needed  on  an  ele¬ 
ment.  Coarsening  leads  to  a  decrease  of  p  by  one. 

A  slight  modification  of  Strategy  B  is  used  to  select  an  initial  grid  for  all  hp- 
strategies.  Beginning  with  an  initial  uniform  grid  of  20  elements  and  uniform  order 
p  =  3,  up  to  ten  iterations  of  the  p-refinement  algorithm  just  discussed  are  executed 
using  Szabo’s  [26]  p-refinement  criterion.  Elements  that  become  linear  during  this 
process  remain  so;  thus,  there  is  no  need  for  extrapolation.  The  iteration  is  continued 
until  either  the  appropriate  tolerance  is  satisfied  on  each  element  or  the  maximum 
number  of  iterations  is  reached,  //-refinement  is  performed  on  the  last  iteration  to  get 
the  initial  grid. 

Initial  conditions  between  meshes  of  different  dimension  are  obtained  by  first 
interpolating  nodal  values  to  the  new  mesh  to  obtain  the  piecewise  linear  approxima¬ 
tion  and  then  using  a  simplified  energy  projection  to  obtain  the  higher-order 
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coefficients.  This  procedure  which  only  requires  local  computations,  and  is  motivated 
by  the  following  lemma. 

Lemma  4.1.  Let  q(x)  £  Hp+l(Cl)  and  Q(x )  £  SAa  with  =p,  i  =  1,  2,  ••• ,N ,  be 
such  that 

Q(xi)  =  q(xi),  (4.3a) 

(q'  -  Q'y')i  =  o,  for  ally  6  sj,(*  (4.3b) 

where 

*/ 

(W',V %  =  J  W'(x)V'(x)dx,  i  =  1,  2,  N.  (4.3c) 

JTi-l 

If  kt  =  xt  -  X:_x  and  k  =  max  x,  -  x._,  then 

IZiZN 

\\q(-)-Q(-)\\i^CkPmilx  (4.3d) 

where  C  is  independent  of  the  grid. 

Proof.  Let  W(x)£  SP/X‘~X,-'\  i.e.,  W (x^)  =  q (x^)  and  W(xi)  =  q(xi).  Setting 
V  =  W  -  Q  £  Sq^X‘  and  using  (4.3b) 

II q'  -  W'Wli  =  II q'  -  Q'Uli  +  II HloV  (4.4) 

Local  Sobolev-space  norms  IH^-  are  defined  as  their  global  counterparts  [23]  but  are 

restricted  to  a  subinterval  (Xf_ltxf).  Neglecting  the  last  term  in  (4.4)  gives, 

\\q'  ~  Q'Woj  ^  II q'  ~  W%,  for  all  W  £  Sp/X'  '  Xi~l\  (4.5) 

Choosing  W (x)  £  SPE  '(x‘  x‘~l)  to  be  be  an  interpolant  and  using  standard  estimates 
[23]  yields 

II q'-Q%&Ck»\\q\\2j  (4.6) 

where  C  is  independent  of  kt.  Conditions  (4.3)  and  use  of  the  Schwarz  inequality 
imply 
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II?  -  Q  llu  S  Ckpto  Wpa.  /  =  1,  2.  W.  (4.7) 

The  result  (4.3d)  follows  from  a  summation  over  i .  □ 

A  strategy  similar  to  that  used  with  BDF  codes  (cf.,  e.g.,  Shampine  and  Gordon 
[24])  is  used  to  eliminate  the  need  for  a  user-prescribed  initial  time  step  and  order. 
The  method  begins  with  an  order  one  SIRK,  and  increases  the  order  by  one  and  dou¬ 
bles  the  time  step  on  successive  steps  until  the  tolerance  is  violated.  Unlike  Shampine 
and  Gordon  [24]  instead  of  ending  this  initial  integration  phase  if  the  temporal  eiror  is 
too  large  on  the  first  time  step,  we  make  a  second  attempt  with  a  time  step  selected  to 
reduce  the  temporal  error  to  1/20  of  the  tolerance.  This  greatly  improved  the  perfor¬ 
mance  of  the  initial  phase  integration  by  allowing  us  to  use  higher-order  methods  in 
time  more  quickly  than  with  Shampine  and  Gordon’s  [24]  algorithm.  To  compare  the 
two  we  also  consider  an  additional  Strategy  ES. 

ES  No  Extra  Startup  Strategy.  Strategy  FP  is  used  with  Shampine  and  Gordon’s  [24] 

initial  integration  strategy. 

//-refinement  is  not  allowed  during  this  initial  phase;  thus,  this  phase  could  end  prema¬ 
turely  if  the  initial  mesh  is  inadequate. 

5.  Computational  Results.  Results  of  four  examples  demonstrate  the  perfor¬ 
mance  of  the  adaptive  h-  and  /^-refinement  procedures.  The  first  three  examples  are 
used  to  compare  the  various  spatial  refinement  strategies;  temporal  enrichment  remains 
invariant  except  for  the  acceptance  or  rejection  of  partial  steps  and  the  initial  integra¬ 
tion  phase.  Experimentation  is  restricted  to  the  /^-refinement  strategy  in  space  (except 
where  noted)  since  the  properties  of  variable  order-variable  step  SIRKs  have  been 
investigated  elsewhere  [8].  The  latter  two  examples  indicate  that  our  general-purpose 
software  can  be  used  to  solve  realistic  problems  with  no  user  intervention.  All  compu¬ 
tations  were  performed  in  double  precision  on  SUN  4/380  and  SPARCstation  II  works¬ 
tations.  All  timings  were  performed  on  SPARCstation  II  workstations. 
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Computational  effort  was  measured  by  either  (i)  the  total  number  of  space-time 
unknowns  including  those  for  the  error  estimates  or  (ii)  the  CPU  time  to  integrate  to 
time  T .  For  Strategy  SX,  the  third  solution  (3.10)  with  k  =  1  was  not  needed  and  this 
is  reflected  in  the  number  of  unknowns  and  the  CPU  time. 

Example  5.1.  Consider  the  linear  heat  conduction  equation 

ut  +f(x,t)  =  uxx,  0<x<l,  t  >  0,  (5.1a) 

with  /  (x  ,t ),  the  initial  data,  and  the  Dirichlet  boundary  conditions  chosen  so  that  the 
exact  solution  is 

u(x,t)  =  2.4tanh75(x  +  1.4/  -  1.4)  -  2tanhl00(x  -  t  -  1/4).  (5.1b) 

Thus,  the  solution  represents  two  steep  wave  fronts  moving  in  opposite  directions  that 
interact  and  pass  through  each  other. 

We  solved  this  problem  for  0  <  t  <.  0.7  with  atol  =  0.05/5*,  k  =  0,  1,  •  ••,  8,  and 
rtol  =  0  using  the  various  h-  and  /zp-refinement  strategies  of  §4.  The  global  errors  in 
the  H 1  norm  at  /  =  0.7  are  presented  as  functions  of  the  number  of  unknowns  in  Fig¬ 
ure  5.1. 

It  is  clear  that  Strategy  FP  is  superior  to  the  //-refinement  Strategies  2,  3,  and  4 
even  for  large  tolerances.  The  //-refinement  Strategy  12  offers  slightly  better  perfor¬ 
mance  as  a  function  of  the  number  of  unknowns  but  it  frequently  failed  to  reduce  the 
error  sufficiently  when  the  tolerance  was  reduced.  Since  //-refinement  only  changes 
the  number  of  elements  by  an  integral  amount,  an  order- 12  method  can  drastically 
reduce  the  local  error;  thus,  satisfying  accuracy  requirements  for  several  tolerances. 
Strategy  p  was  not  able  to  solve  (5.1)  when  atol  =  0.05/57  and  0.05/58. 

The  adaptive  /zp-refinement  strategies,  shown  at  the  right  and  bottom  of  Figure 
5.1,  all  produce  comparable  results  for  the  larger  tolerances  and  errors.  Lower-order 
SIRKs  are  used  with  larger  tolerances  and  any  partial  steps  that  are  taken  are  only  a 
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Figure  5.1.  Global  errors  in  the  Z/1  norm  at  /  =0.7  for  Example  5.1  as  a 
function  of  the  number  of  unknowns  for  the  Strategies  FP,  2,  3,  4,  12,  p, 
(upper  left),  and  B,  FP,  SX,  NP,  NG,  and  ES  (upper  right).  Global  errors  in 
the  H  norm  at  t  =  0.7  for  Example  5.1  as  a  function  of  the  CPU  time  for 
the  Strategies  B,  FP,  SX,  NP,  NG,  and  ES  (bottom). 


Figure  5  2.  Grid  and  order  used  in  solving  Example  5.1  with  Strategy  FP 
with  a  to  l  -  0.01.  Temporal  order  is  indicated  by  greyscale,  and  spatial  ord- 
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1,0 
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0 

2,0 

1,0 

1,0 

1,0 

1,0 
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69,1 

0.01/56 

2,0 

1,0 

1,0 

1,0 

1,0 

■ 

158,20 

7,3 

0.01/57 

2,0 

1,0 

1,0 

1,0 

1,0 

mWm 

188,16 

168,12 

o 

o 

ui 

00 

2,0 

1,0 

1,0 

1,0 

1,0 

10,0 

114,1 

361,21 

Table  5.1.  The  number  of  time  steps  and  the  number  of  partial  time  steps 
taken  as  a  function  of  the  number  of  stages  and  tolerance  when  Strategy  FP 
is  applied  to  Example  5.1. 


atol 

Strategy 

g||||j 

m 

l 

llll^ 

0.01/55 

FP 

0.7% 

0.6% 

0.4% 

0.6% 

0.3% 

0.4% 

0.2% 

0.6% 

NG 

9.4% 

6.2% 

5.3% 

3.0% 

3.6% 

4.9% 

3.3% 

2.1% 

Table  5.2.  Percentage  of  rejected  steps  when  Strategies  FP  and  NG  are  ap¬ 
plied  to  Example  5.1. 

small  fraction  of  a  potential  full  step.  However,  when  the  higher-order  SIRKs  are  used 
with  the  smaller  tolerances,  there  is  a  clear  advantage  of  the  partial-step  strategies 
(e.g.,  FP  and  SX)  relative  to  Strategy  NP.  Partial  steps  with  the  higher-order  SIRKs 
may  be  as  much  as  40%  of  a  potential  full  step.  Table  5.1  indicates  the  number  of 
time  steps  and  the  number  of  partial  time  steps  that  were  taken  at  each  SIRK  order  and 
each  tolerance  for  Strategy  FP.  These  results  verify  that  the  higher-order  SIRKs  are 
taking  a  sizable  number  of  partial  time  steps. 


At  the  smaller  tolerances,  Strategies  FP  and  SX  are  clearly  superior  to  B;  thus, 
the  ability  to  predict  a  future  mesh  is  a  definite  advantage.  Strategy  ES  is  slightly  less 
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efficient  than  FP  at  the  smallest  tolerances,  but  is  otherwise  comparable.  Strategy  NG 
has  an  undesirable  erratic  behavior  as  the  tolerance  is  varied.  Finally,  the  two  com¬ 
plexity  measures  are  comparable  with  Strategy  SX  being  the  most  efficient  relative  to 
CPU  time. 

The  percentage  of  rejected  steps  is  shown  as  a  function  of  tolerance  for  Strategies 
FP  and  NG  in  Table  5.2.  Once  again,  the  ability  to  predict  future  grids  is  valuable  and 
has  reduced  the  number  of  rejected  steps. 

The  grid  used  to  compute  the  solution  using  Strategy  FP  with  atol  =  0.01  is 
shown  in  Figure  5.2.  The  greyscale  indicates  the  temporal  order  used  and  the  height 
( P  axis)  indicates  the  spatial  order.  The  grid  and  method  order  both  track  the  steep 
fronts.  After  interaction,  the  fine  grids  separate  following  the  individual  waves.  In 
this  example  the  temporal  order  quickly  climbs  to  five  (during  the  initial  phase)  and 
remains  there  throughout  the  integration. 

Example  5.2.  Consider  Burgers’  equation 

ut  +  uux  =  su^,  0  <  x  <  1,  t  >  0,  (5.2a) 

where  e  and  the  initial  and  Dirichlet  boundary  conditions  are  chosen  so  the  exact  solu¬ 
tion  is 

u(x,t)  -  l-2tanh81(x  -  t  +  1.3).  (5.2b) 

This  solution  is,  once  again,  a  steep  wave  traveling  in  the  positive  x  direction  as  time 
increases. 

We  solved  this  problem  for  0  <  t  <,  1.0  with  atol  =  5x  10-*,  k  =  0,  1,  7,  and 

rtol  =  0  using  Strategy  B  and  the  time  extrapolation  strategies  described  in  §4.1.  The 
error  in  the  Hl  norm  is  presented  as  a  function  of  the  number  of  unknowns  in  Figure 
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5.3.  The  three  Strategies  BBC,  MF  and  CO  are  comparable  for  all  but  the  smallest 
tolerance  where  BBC  performs  slightly  better.  Strategies  IC  and  MB  performed 
significantly  worse  than  the  other  three.  The  optimal  Newton  initial  guess  strategy  is 
dependent  on  both  the  problem  and  on  the  adaptive  method  used  to  solve  it. 

We  also  solved  (5.2)  for  the  same  time  interval  and  tolerances  with  the  hp- 
refinement  strategies  of  §4  and  the  “best”  temporal  extrapolation  strategy.  In  the  case 
of  Strategies  B,  FP,  SX,  and  NG,  the  best  Newton  guess  strategy  is  BBC  while  NP 
works  best  with  MF.  The  results,  given  in  Figure  5.4,  indicate  that  Strategies  FP  and 
SX  are  again  clearly  superior  to  NP  and  B.  Strategy  NG  uses  fewer  unknowns  for 
small  tolerances  but,  the  difference  between  it  and  Strategy  FP  as  measured  in  CPU 
time  is  not  significant.  As  in  Example  5.1,  Strategy  NP  is  comparable  to  FP  and  SX 
for  the  large  tolerances  but  becomes  less  efficient  for  small  tolerances.  Thus,  accept¬ 
ing  partial  steps  seems  to  be  beneficial  when  solving  partial  differential  equations  to 
high  accuracy.  Strategies  ES  and  FP  produce  identical  results  for  this  problem.  In 
Table  5.3,  we  again  verify  that  Strategy  NG  has  rejected  a  larger  percentage  of  steps 
than  Strategy  FP. 

Example  5.3.  Consider  the  Brusselator  problem  with  diffusion  [15] 

ut  -  1  -  u2v  +  4.4m  =  eiixx,  (5.3a) 

vt  -  3.4m  +  u2v  =  ev^,  0  <  x  <  1,  t  >  0,  (5.3b) 

m(x,0)  =  0.5,  v(x,0)  =  1  +  5x  -  tanh(20x)/4  -  tanh(20(x  -  l))/4,  0  <  x  <  1, 

(5.3c) 

ux(0,t)  =  ux(l,t)  =  vx(0,0  =  vjC(l,/)  =  0,  t>  0.  (5.3d) 

The  boundary  data  prescribed  for  v  (x  ,0)  is  consistent  with  (5.3d). 
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Figure  5.3.  Global  errors  in  the  Z/1  norm  at  t  -  1.0  as  a  function  of  the 
number  of  unknowns  when  the  Newton  initial  guess  Strategies  BBC,  MF, 
CO,  IC,  and  MB  are  applied  to  Example  5.2. 
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Figure  5.4.  Global  errors  in  the  Z/1  norm  at  t  =  1.0  as  a  function  of  the 
number  of  unknowns  (left)  and  CPU  time  (right)  when  Strategies  B,  FP  SX 
NP,  and  NG  are  applied  to  Example  5.2. 
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atol 
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0.01 
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ul 

0.01/52 

0.01/53 

0.0 1/54 

0.01/55 

0.01/56 

0.01/57 

0.01/58 

FP 

10.3% 

7.8% 

5.5% 

25% 

3.5% 

0.7% 

1.2% 

2.4% 

0.5% 

NG 

11.6% 

10.6% 

10.0% 

6.0% 

4.9% 

4.3% 

6.0% 

5.1% 

2.4% 

Table  5.3.  Percentage  of  rejected  steps  when  Strategies  FP  and  NG  are  ap¬ 
plied  to  Example  5.2. 

We  solved  this  problem  for  0  <  /  £  12.6  with  e  =  0.002  and 
atol1  =  rtol1  -  0.05/5*,  i  =  1,  2,  k  =  0,  1,  ",  8,  using  Strategies  B,  FP,  NG,  SX,  and 
ES  with  Strategy  BBC.  As  an  “exact  solution,”  we  used  the  results  obtained  from 
Strategy  FP  with  atol1  =  rtol1  =  0.05/513,  i  =  1,  2.  The  difference  between  this  exact 
solution  and  a  Strategy  FP  solution  with  atol1  =  rtol1  =  0.05/5 12,  i  =  1,2,  is 
5.38  x  10“9  in  Hl.  The  H1  error  is  shown  as  a  function  of  the  number  of  unknowns  in 
Figure  5.5.  The  solution  components  at  t  =  12.6  exhibit  steep  fronts  as  shown  in  Fig¬ 
ure  5.6.  Unlike  the  previous  examples.  Strategy  B  performed  better  than  Strategy  FP 
for  the  larger  tolerances;  however,  it  exhibits  some  erratic  behavior  as  the  tolerance  is 
varied.  Strategy  SX  also  exhibits  some  erratic  behavior  but  is  clearly  superior  to  Stra¬ 
tegy  FP  for  the  smaller  tolerances.  Strategies  FP  and  SX  are  superior  to  and  less 
erratic  than  Strategies  B,  NG,  and  NP  at  the  smaller  tolerances.  Apparently,  using 
advanced  grid  information  to  determine  /z-refinement  has  a  smoothing  effect  on  the 
convergence  rate  of  the  method.  Strategies  FP  and  SX  also  benefit  from  our  modified 
initial  integration  phase. 

The  grid  and  method  order  used  to  calculate  the  solution  on  0  <  t  <,  12.6  of  (5.3) 
using  Strategy  FP  with  BBC  and  atol1  =  rtol1  =  0.01,  i  =  1,  2,  is  shown  in  Figure  5.7. 
As  in  Example  5.1,  the  finest  grids  and  highest  order  elements  are  concentrated  near 
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Figure  5.7.  Grid  and  order  used  for  Example  5.3  using  Strategy  FP  with 
BBC  and  atol1  =  rtol1  =0.01,  i  =1,2.  Temporal  order  is  indicated  by 
greyscale  and  spatial  order  by  height. 
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the  steep  fronts. 

Example  5.4.  Consider  the  model  of  shear  band  formation  introduced  by  Drew 


and  Flaherty  [12] 

ut  -  v  =  0,  v,  =  [G(T)ux]x  +  v^/Re,  (5.4a,b) 

Tt  -  ( vx)2IRe  =  TjaKPrRe),  0  <  x  <  1,  /  >  0,  (5.4c) 

G(T)  =  1/2[(1  +  Goo)  -  (1  -  G oo) tanh((7 -Tm )/AT )] ,  (5.4d) 

u(x,0)  =  v(x,0)  =  T(x,0)  =  0,  0  < x  <  1,  (5.4e) 

v  (0,t )  =  0,  v (1, t)  =  V(t),T (0,0  =  0,  7(1,0  =  0,  r  >  0.  (5.4f) 


This  system  describes  the  simple  shearing  of  a  slab  of  material  having  unit  thickness. 
The  variables  u ,  v ,  and  T  denote  the  displacement,  velocity,  and  temperature  of  the 
material  point  at  the  position  x  and  time  t .  The  shear  modulus  G  depends  on  tem¬ 
perature  according  to  (5.3d);  thus,  when  0  <  AT  1  the  material  undergoes  a  phase 
transition  at  the  temperature  Tm  with  the  shear  modulus  abruptly  changing  from  unity 
to  Goo.  One  edge  (x  =  0)  of  the  slab  is  held  fixed  while  the  other  (x  =  1)  is  subjected 
to  a  shearing  velocity  V(t).  When  F(r)  varies  slowly,  the  velocity  v  is  approximately 
a  linear  function  of  x;  however,  if  V  (t)  varies  rapidly  and  the  Prandtl  number  Pr  and 
Reynolds  number  Re  have  appropriate  values,  the  deformation  will  be  localized  in  a 
narrow  region  called  a  shear  band.  Drew  and  Flaherty’s  [12]  model  (5.4)  omits  physi¬ 
cal  considerations  which  have  been  included  in  studies  by  Wright  and  Walter  [28]  and 
Walter  [27];  however,  (5.4)  serves  to  illustrate  that  difficult  nonlinear  systems  can  be 
solved  without  a  priori  knowledge  of  solution  behavior. 


Consider  a  problem  [12]  with 
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V(t)=  1 


V0(tlr), 

V0 , 

W  "  O/r, 

0, 


0  <  t  <  r 
r  <  d  -  r 
d  -  r  <  t  <  d 
d  <  t. 


(5.4g) 


The  initial  shearing  velocity  initiates  a  wave  that  propagates  through  the  solid  medium. 
The  unloading  at  i  =  d  -  r  sends  a  second  wave  through  the  material.  The  passage  of 
these  waves  generates  heat  which  reduces  the  shear  modulus.  Sufficient  heat  can 
cause  a  phase  transition  to  occur  which  can  localize  the  deformation.  In  order  to  illus¬ 
trate  this,  we  performed  a  computation  with  Re  =  100,  Pr  =  50,  GM  =  0.05, 
Tm  =  0.03,  AT  =  0.01,  V0  =  0.5,  d  =  1.5,  and  r  =  0.05  and  solved  (5.4)  for 
0  ^  t  <;  3.24  with  tolerances  atol1  =  0.01  and  rtol1  =  0,  i  =  1,  2,  3,  using  Strategy  FP 
with  BBC.  Solutions  are  shown  as  functions  of  *  for  t  -  0.18,  0.47,  0.78,  1.16,  1.46, 
1.69,  1.95,  2.18,  2.65,  and  3.08  in  Figure  5.8.  The  space-time  grid  used  for  this  calcu¬ 
lation  is  shown  in  Figure  5.9.  With  this  data,  the  temperature  rises  rapidly  and  a  shear 
band  forms  in  a  layer  adjacent  to  the  loaded  edge  ( x  =  1).  The  grid  is  concentrated  in 
the  shear  band  region  and  follows  its  evolution  into  the  domain.  Initial  orders  increase 
quickly  and  then  decrease  where  they  are  not  needed  for  accuracy. 


6.  Discussion.  We  present  and  compare  several  adaptive  /ip-refinement  strategies 
for  solving  parabolic  systems.  This  investigation  is  much  more  extensive  than  previ¬ 
ous  studies  [3,  19]  which  used  larger  tolerances  and  no  grid  prediction.  It  also  unifies 
spatial  and  temporal  enrichment  to  a  much  higher  degree  than  was  done  previously. 
Spatial  Zip -refinement  provides  significant  improvements  in  efficiency  relative  to  h- 
refinement  with  either  low-  or  moderate- order  methods. 


The  differences  in  efficiency  between  the  various  ftp-refinement  strategies  in 
space  are  not  as  great.  Taking  partial  steps  offers  some  advantages  with  small  toler- 


-  31  - 


Figure  5.8.  Displacement  (upper  left),  velocity  (upper  right),  and  tempera¬ 
ture  (bottom)  as  functions  of  x  for  Example  5.4  at  several  times  using  FP 
with  BBC  and  atol1  -  0.01,  rtoV  =  0,  i  =  1,  2,  3.  6 


ances  and  high-order  methods.  Grid  prediction  (Strategies  B,  FP,  and  SX)  leads  to  a 
smoother  relationship  between  the  error  and  tolerance  than  no  prediction  (Strategy 
NG).  Although  Strategy  SX  outperformed  Strategy  FP  in  almost  every  case,  indicating 
that  extrapolation  of  low-order  solutions  can  produce  adequate  indications  of  high- 
order  error,  it  did  exhibit  some  erratic  behavior  for  large  tolerances  in  Example  5.3. 
Further  investigation  is  needed  to  clarify  which  is  the  better  approach.  In  most  cases 
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the  number  of  unknowns  and  CPU  time  were  comparable  measures  of  work.  How¬ 
ever,  there  were  exceptions  in  Example  5.1  with  Strategy  SX  and  in  Example  5.2  with 
Strategy  NP. 

The  three  strategies  for  computing  initial  guesses  for  Newton’s  method,  BBC, 
MF,  and  CO,  were  comparable  with  the  best  Newton  guess  strategy  being  dependent 
upon  the  adaptive  strategy  (cf.  Example  5.2)  and,  clearly,  the  problem  as  well.  The 
grids  and  orders  selected  according  to  Strategy  FP  successfully  track  important  solution 
features.  Examples  5.3  and  5.4  indicate  that  the  adaptive  /zp-refinement  strategies  can 
robustly  solve  some  difficult  nonlinear  systems  with  no  user  intervention. 

Although  no  attempt  was  made  to  control  the  global  temporal  error,  it  appears 
that  we  have  done  so  to  a  high  degree.  With  stiff  ordinary  differential  systems,  such 
as  those  arising  here,  local  errors  do  not  accumulate  [13,  17];  hence  control  of  global 
errors  is  accomplished  simply  by  controlling  local  errors. 

Several  issues  remain.  An  alternative  spatial  error  estimation  strategy  which 
involves  solving  local  elliptic  problems  rather  than  local  parabolic  problems  has  been 
implemented  in  a  MOL  code  which  uses  BDF  formulas  in  time  [3].  Reported  gains  in 
efficiency  could  be  even  greater  with  SIRKs  since  spatial  error  estimates  are  needed 
only  at  certain  stages  for  the  particular  adaptive  strategy  used.  Theoretical  analysis  of 
these  estimates  with  SIRKs  was  done  by  Moore  in  [21]  to  show  that  they  produce 
asymptotically  correct  results.  A  comparison  of  SIRK  and  BDF  methods  must  be 
done.  The  BDF-based  methods  may  outperform  SIRKs  at  the  larger  tolerances  but  this 
may  not  be  so  for  smaller  tolerances.  The  strategies  and  software  developed  herein 
will  be  most  helpful  in  developing  multi-dimensional  techniques. 


Figure  5.9.  Grid  on  0  <;  t  <.  3.24  (left)  and  order  on  0  ^  t  <.  0.38  (right)  used 
for  Example  5.4  using  Strategy  FP  with  BBC  and  atol1  =  0.01,  rtol1  =  0, 
i  =  1,  2,  3.  Temporal  order  is  indicated  by  greyscale  and  spatial  order  bv 
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